%% Figure 2A: Color plots for 3 devices. 
% D0A1393 
a = 0;
b = 15;

figure('position', [10 10 (b-a)*21/0.8 400]);
set(gca, 'position', [0.1,0.1,0.8,0.8]);
hold off;
plot(D0A1393.Ic(:, 2)*D0A1393.RnTc2 * 1000 , D0A1393.Ic(:, 1), 'color', [0,0,0]);
xlim([-b, b]);
ylim([10, 90]);
xlabel('I R (mV)');
ylabel('Temperature (K)');
%title('Interface 2; 0 degrees');

figure('position', [10 10 (b-a)*21/0.8 400]);
hold off;
set(gca, 'position', [0.1,0.1,0.8,0.8]);
sdVdI = zeros(size(D0A1393.CdVdI));
for i = 1:53
    sdVdI(i, :) = smooth(D0A1393.CdVdI(i, :), 3);
end
pcolor(D0A1393.I(:, 2:end)*D0A1393.RnTc2 * 1000, D0A1393.TSample(:, 2:end), sdVdI / D0A1393.RnTc);
shading interp;
caxis([0, 10]);
map = GenerateColorScaleJunctions();
colormap(map);
hold all;


xlim([-b, b]);
ylim([10, 90]);
 
xlabel('I R (mV)');
ylabel('Temperature (K)');
title('Interface 2; 0 degrees');

% D31A417
a = -5.4;
b = 5.4;
figure('position', [10 10 (b-a)*10.5/0.8 400]);
set(gca, 'position', [0.1,0.1,0.8,0.8]);
hold off;

pcolor(D31A417.I(:, 2:end) * D31A417.RnTc2 * 1000, D31A417.TSample(:, 2:end), D31A417.CdVdI / D31A417.RnTc2);
shading interp;

xlim([a, b]);
ylim([10, 90]);
caxis([0, 10]);
map = GenerateColorScaleJunctions();
colormap(map);

xlabel('I R (mV)');
ylabel('Temperature (K)');
title('Interface 15A; 27 degrees');

% D449A123 - Interface 19 B
a = -2
b = 2;
figure('position', [100 100 (b-a)*20/0.8*2 400]);
set(gca, 'position', [0.1,0.1,0.8,0.8]);
hold off;
pcolor(D449A123.I(:, 1:end) / D449A123.Area * 1e5, D449A123.TSample(:, 1:end), D449A123.dV(:, 1:end) ./ D449A123.dI(:, 1:end) / D449A123.RnTc2);

shading interp;

ylim([10, 90]);
caxis([0, 10]);
map = GenerateColorScaleJunctions();
colormap(map);

ylabel('Temperature (K)');
title('Interface 7 B; 43 degrees');

%% Figure 2B: IcRn vs Angle at 30K
figure(1);
clf;
x = 0:90;
i = 7;

errorbar(Ic2RnMasterTable([2, 5, 7:27], 1), Ic2RnMasterTable([2, 5, 7:27] , i) * 1000, nanstd(Ic2RnMasterTable([2, 5, 7:27] , i-3:i+3), 0, 2) * 1000, 's', 'MarkerFaceColor', [39,142,255] / 255, 'MarkerEdgeColor', [39,142,255] / 255, 'color', [39,142,255] / 255 ); 
xlim([0, 90]); 
hold all;
errorbar(mod(Ic2RnMasterTable([3, 4, 6], 1), 90), Ic2RnMasterTable([3, 4, 6], i) * 1000, nanstd(Ic2RnMasterTable([3, 4, 6], i-3:i+3), 0, 2) * 1000, 's',  'MarkerEdgeColor', [39,142,255] / 255, 'color',  [39,142,255] / 255);

ylabel('IcRn (mV)');
ylim([0, 25]);

colors = [0,0,0; 0,0,255; 39, 142, 255; 22,216,105] / 255;

i = 25;
errorbar(Ic2RnMasterTable([2, 5, 7:27], 1), Ic2RnMasterTable([2, 5, 7:27] , i) * 1000, nanstd(Ic2RnMasterTable([2, 5, 7:27], i-3:i+3), 0, 2) * 1000, 'o', 'MarkerFaceColor', [237,62,125]/255, 'MarkerEdgeColor', [237,62,125]/255, 'color', [237,62,125]/255); 
xlim([0, 90]); 
set(gca, 'XTick', [0, 22.5, 45, 67.5, 90]); 
hold all;
errorbar(mod(Ic2RnMasterTable([3, 4, 6], 1), 90), Ic2RnMasterTable([3, 4, 6], i) * 1000, nanstd(Ic2RnMasterTable([3, 4, 6], i-3:i+3), 0, 2) * 1000, 'o',  'MarkerEdgeColor', [237,62,125]/255, 'color', [237,62,125]/255);
plot(x, 17.22*abs(cos(2*x*pi/180)), 'color', 0.5*[1,1,1])
ylabel('IcRn (mV)');
ylim([0, 25]);
pause(1);

x = 0:90;

colors = [0,0,0; 0,0,255; 39, 142, 255; 22,216,105] / 255;

%% Figure 2C Angle-Dependent Ic Suppression at Low T
% Plot Ic*Rn(Tc) vs T for subset of devices in one figure.
varList = whos('D*');

%Take only 0, 27, 31, 39, 43 degrees data.
%varList = varList([3, 11, 12, 13, 15, 9]);
%Angles =         [0, 14, 119, 39, 43]
colors = [8, 73, 255; 
    65, 174, 229;
    65, 229, 174;
	255, 85, 166;
132, 68, 224;
206, 0, 84;
255,64,210;
102,46,255] / 255

figure(2);
hold off;
title('Critical Current Density over All Devices');

%colors = jet(length(varList));
% %Agl  [0, 14, 119, 39, 43]
% varList = varList([3,  8,   5, 13, 15]);
% Ta = [4 , 5 , 4,  4,   4];
% Tb = [25, 18, 26, 29,  35];

%Agl  [14, 119, 39, 43]
varList = varList([1,   3, 11, 20]);
Ta = [5 , 4,  4,   4];
Tb = [18, 26, 29,  35];

%Agl  0 ,  0, 84, 14, 76, 15, 70,157, 27,119, 31, 57,129, 39,  39,  43, 134
%Ta = [7 , 5 , 10,  5,  4,  6,  4,  4,  4,  7, 10,  5,  4,  4,   6,  28,  46]
%Tb = [25, 25, 25, 18, 12, 30, 10, 16, 13, 26, 24, 13, 29, 29,  24,  52,  56]

for i = [1,2,3]
    
    varList(i).name
    eval(strcat('data = ', varList(i).name, ';'));
    plot(data.Ic2(:, 1), data.Ic2(:, 2) * data.RnTc2 * 1000, 'o', 'markersize', 4, 'color', colors(i, :), 'linewidth', 1);
    hold all
    
    [~, idxA] = min(abs(data.Ic2(:, 1) - Ta(i)));
    [~, idxB] = min(abs(data.Ic2(:, 1) - Tb(i)));
    if idxA > idxB
        ans = idxB;
        idxB = idxA;
        idxA = ans;
    end
    [m, b, r2] = FitLine(data.Ic2(idxA:idxB, 1), data.Ic2(idxA:idxB, 2) * data.RnTc2 * 1000);
          
    plot(Ta(i):Tb(i) + 7, m(1) * (Ta(i):Tb(i)+7) + b(1), 'color', colors(i, :), 'linewidth', 1.5);
       
    xlabel('Temperature (K)');
    ylabel('Ic Rn (mV)');

end

    e = 1.602176e-19;
    IcRnD = (pi * 0.016595 * Gap.dWave / (2)) .* tanh(0.016595*Gap.dWave ./ (2*8.6173e-5 * 90 * Gap.x));

    plot(90 * Gap.x, 1000 * 0.65 * IcRnD * 0.6 );


i = 4
    varList(i).name
    eval(strcat('data = ', varList(i).name, ';'));
    plot(data.Ic2(:, 1), data.Ic2(:, 2) * data.RnTc2 * 1000 * 5, 'o', 'markersize', 4, 'color', colors(i, :), 'linewidth', 1.5);
    
    [~, idxA] = min(abs(data.Ic2(:, 1) - Ta(i)));
    [~, idxB] = min(abs(data.Ic2(:, 1) - Tb(i)));

    if idxA > idxB
        ans = idxB;
        idxB = idxA;
        idxA = ans;
    end
    [m, b, r2] = FitLine(data.Ic2(idxA:idxB, 1), data.Ic2(idxA:idxB, 2) * data.RnTc2 * 1000);

    plot(Ta(i):Tb(i) + 5, 5*(m(1) * (Ta(i):Tb(i)+5) + b(1)), 'color', colors(i, :), 'linewidth', 1.5);


xlim([0, 90]);
ylim([0, 16]);

clear ans;
clear data;
clear i;
clear list;
clear varList;


%% Figure 2D, 2E, S12, 
varList = whos('D*');

%Order runs by distance from 45*. Use abs(45 - mod(angle, 90)).
idx = [2	1	26	5	25	24	6	7	8	3	9	23	4	11	12	10	13	17	14	15	22	16	19	18	21	20];
varList = varList(idx);



figure(1);
clf;
figure(2);
clf;
figure(3);
clf;
figure(4);
clf;
figure(5);
clf;
figure(6);
clf;

T = 1:90;

Ic2Matrix = NaN * ones(length(T), 3*length(varList));
angles = zeros(3*length(varList), 1);


%Order by \theta folded to [0, 45]
%Agl  0,   0, 84, 14, 76, 70, 157, 23, 27, 119, 31, 57, 129, 39, 39, 39.5, 40.3, 43.0, 43, 43.2, 46.3, 43.8, 44.6,44.6,45.2, 44.9
%a    0,   0,  6, 14, 14, 20, 23,  23, 27,  29, 31, 33,  39, 39, 39, 39.5, 40.3, 43.0, 43,       43.7, 43.8,           44.8, 44.9
Ta = [4 ,  5, 10,  5,  4,  4,  4,   4,  4,   5,  5,  4,   4,  4,  4,    4,    4,    4,  4,  4,      4,   4,     4,   4,   4,    4];%range of fit for Tm
Tb = [25, 25, 25, 18, 30, 10, 13,  35, 12,  26, 24, 15,  29, 32, 22,   35,   35,   35, 35, 35,     35,  35,    35,  35,  30,   35];
off =[-5,  3,  10, 20, 27, 31, 37,  45, 50, 52, 54, 63,  65, 69, 73,   74,   77,   75, 78, 82,     85,  86,    88,  90,  92,   95];%Fig S12 trace offset
fac = [1,  1,  1,  1,  1,  1,  1,   1,  1,   1,  1,  1,   1,  1,  1,    1,    1,    1,  1,  1,      1,   1,     1,   1,   1,    1];
Tm = [ 0,  0,  0,  0,  0,  0,  0,   0, 12,  26, 26, 15,  32, 30, 24,    0,    0,   0, 58,  0,      0,   0,     0,   0,   0,    0]; %Start point of fit for Tm.


colors = jet(length(varList));
fitColor = [1,0,0];

for i = 1:length(Ta)
    eval(strcat('data = ', varList(i).name, ';'));
	disp([num2str(i), ' ', varList(i).name]);
        
    %plot the data first.
    figure(1);
    plot(data.Ic2(:, 1), fac(i) * data.Ic2(:, 2) * data.RnTc2 * 1000 - off(i), 'o', 'color', colors(i, :), 'markersize', 4);
    hold all;
    plot([85, 90], [0,0] - off(i), 'color', colors(i, :), 'linewidth', 2)

    %now fit the data in the T range selected.
    [~, idxA] = min(abs(data.Ic2(:, 1) - Ta(i)));
    [~, idxB] = min(abs(data.Ic2(:, 1) - Tb(i)));
    if idxA > idxB
        ans = idxB;
        idxB = idxA;
        idxA = ans;
    end
    [m, b, r2] = FitLine(data.Ic2(idxA:idxB, 1), data.Ic2(idxA:idxB, 2) * data.RnTc2 * 1000);
    
    %overlay the fit on the data.
    figure(1);
    plot(data.Ic2(idxA:idxB, 1), fac(i) * data.Ic2(idxA:idxB, 2) * data.RnTc2 * 1000 - off(i), 'o', 'markersize', 4, 'color', colors(i, :), 'MarkerFaceColor', colors(i, :));
    lineRange = (min(data.Ic2(idxA:idxB, 1)) - 3):max((data.Ic2(idxA:idxB, 1)) + 3);
    
    plot(lineRange, fac(i) * (m(1) * lineRange + b(1)) - off(i),  'color', colors(i, :), 'linewidth', 2);
    
    %plot fitted slope.    
    figure(2);
    if data.Angle < 90
        plot(data.Angle, (m(1)), '^', 'color', [237,62,125] / 255, 'MarkerFaceColor', [237,62,125] / 255);
    else
        plot(data.Angle - 90, (m(1)), '^', 'color', [237,62,125] / 255);
    end
    hold all;
        
    %plot slope vs. reflected angle in range 0 - 45 degrees.
    figure(3);


        plot(mod(data.Angle, 90), (m(1)), '^', 'color', [237,62,125] / 255, 'MarkerFaceColor', [237,62,125] / 255);
 
    hold all;
    %pause(2);
    
    figure(4);
    if data.Angle <= 45
        errorbar(-abs(45 - mod(data.Angle, 90)), m(1), m(2), 'd', 'color', [65,174,229] / 255, 'MarkerFaceColor', [65,174,229] / 255, 'MarkerSize', 8);
    elseif data.Angle < 90
        errorbar(-abs(45 - mod(data.Angle, 90)), m(1), m(2), 'o', 'color', [63,232,163] / 255, 'MarkerFaceColor', [63,232,163] / 255, 'MarkerSize', 8);
    elseif data.Angle < 135
        errorbar(-abs(45 - mod(data.Angle, 90)), m(1), m(2), 'd', 'color', [65,174,229] / 255, 'MarkerSize', 8);
    elseif data.Angle < 180
        errorbar(-abs(45 - mod(data.Angle, 90)), m(1), m(2), 'o', 'color', [63,232,163] / 255, 'MarkerSize', 8);
    end
    
    hold all;
    
    %fit Ic2Rn near T_M with a Gaussian, to see where it is and how wide it is. 
    if Tm(i) ~= 0
        [~, idx1] = min(abs(Tm(i) - data.Ic2(:, 1)-15 * fac(i))); 
        [~, idx2] = min(abs(Tm(i) - data.Ic2(:, 1)+15*fac(i)));
        if idx1 < idx2
            idx = idx1:idx2;
        else
            idx = idx2:idx1;
        end
        [A, mu, sigma, m, b] = FitPeak(data.Ic2(idx, 1), 1000 * data.Ic2(idx, 2) * data.RnTc2);
        
        Tm(i) = mu(1);
        TmWidth(i) = sigma(1);
        
        figure(1);
        plot(data.Ic2(idx, 1), fac(i) * (A / (sigma * sqrt(2*pi)) * exp(-1/2 * (data.Ic2(idx, 1)-mu).^2 / sigma.^2) + m*data.Ic2(idx, 1) + b) - off(i));
    else
        %set to zero.
        Tm(i) = 0;
        TmWidth(i) = 0;
    end
    
    figure(5);
    
    if (data.Angle < 90)
        errorbar(data.phi, Tm(i), TmWidth(i), 'd', 'MarkerSize', 10, 'color',  [237,62,125] / 255, 'CapSize', 0, 'MarkerFaceColor',[237,62,125] / 255 );
    else
        errorbar(data.phi, Tm(i), TmWidth(i), 'd', 'MarkerSize', 10, 'color',  [237,62,125] / 255, 'CapSize', 0);
    end
     hold all;
    
    figure(6);
        
    if (data.Angle < 45)
        errorbar(-abs(45 - mod(data.Angle, 90)), Tm(i), TmWidth(i), 'd', 'MarkerSize', 8, 'color',  [65,174,229] / 255, 'CapSize', 0, 'MarkerFaceColor',[65,174,229] / 255 );
    elseif (data.Angle <= 90)
        errorbar(-abs(45 - mod(data.Angle, 90)), Tm(i), TmWidth(i), 'o', 'MarkerSize', 8, 'color',  [63,232,163] / 255, 'CapSize', 0, 'MarkerFaceColor',[63,232,163] / 255 );
    elseif (data.Angle <= 135)
        errorbar(-abs(45 - mod(data.Angle, 90)), Tm(i), TmWidth(i), 'd', 'MarkerSize', 8, 'color',  [65,174,229] / 255, 'CapSize', 0);
    else
        errorbar(-abs(45 - mod(data.Angle, 90)), Tm(i), TmWidth(i), 'o', 'MarkerSize', 8, 'color',  [63,232,163] / 255, 'CapSize', 0);
    end
     hold all;
    
end

figure(1);
xlabel('Temperature (K)');
ylabel('IcRn (mV)');
xlim([0, 90]);
ylim([-97, 25]);

figure(2);
xlabel('\phi');
ylabel('d[IcRn]/dT');
plot([0, 90], [0, 0]);
xticks([0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2] * 90);

figure(3);
xlabel('\alpha');
ylabel('|d[IcRn]/dT|');
plot([0, 90], [0,0]);
xlim([0, 90]);
xticks([0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2] * 90);


figure(4);
xlabel('| 45 - \phi |');
ylabel('d[IcRn]/dT (mV/K)');
plot([-45, 0], [0,0]);
xlim([-45, 0]);
xticks([-1, -0.75, -0.5, -0.25, 0] * 45);

figure(5);
xlabel('\phi');
ylabel('Tm (K)');
xlim([0, 90]);
xticks([0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2] * 90);

figure(6);
xlabel('|45 - \phi|');
ylabel('Tm (K)');
xticks([-1, -0.75, -0.5, -0.25, 0] * 45);







%% S6: IV Shape Overview.
figure(20); clf; figure(30); clf;
%D119A736 has best data to high voltage.
data = D119A736;
midpt = ceil(length(data.V)/2);

ft = fittype('CalcCurrent(x, Rn, Delta, T)', 'problem', 'T', 'independent', {'x'}, 'coefficients', {'Rn', 'Delta'});

%specify list of traces to plot here: These are row indices in each data
%file. Each row is a scan at a different temperature.
list = [5, 10, 20, 40, 60, 80];   

colors = [21, 89, 209;
        221, 199, 255;          
		  242, 22, 112];
      


colors = interp1([0; 0.5; 1], colors, linspace(0, 1, length(list))) / 255;

for i = 1:length(list)
	[~, idx] = min(abs(list(i) - data.Ic2(:, 1)));	

	figure(20);
    	plot(-(data.V(idx, :) - data.V(idx, midpt)) * 1000, -data.I(idx, :) / data.Area / 1000 * 1e8,  'color', colors(i, :), 'linewidth', 2, 'DisplayName', strcat(num2str(round(mean(data.TSample(idx, :))), 2), 'K') );
		hold all;
	
   	figure(30);
		plot(-(data.V(idx, :) - data.V(idx, midpt)) * 1000, (data.dV(idx, :) ./ data.dI(idx, :)), '.',  'color', colors(i, :), 'linewidth', 2, 'DisplayName', strcat(num2str(round(mean(data.TSample(idx, :))), 2), 'K') );		hold all;
	
end

% fit the low T data to theory.
i = 83;
figure(10);
%plot I / area in [uA / um^2], and V (and correcting for instrument offset) in [mV]
semilogy((-(data.V(i, :) - data.V(i, midpt))) * 1000, -data.I(i, :) / data.Area / 1000 * 1e8, '.', 'color', colors(1, :))

[~, endIdx] = min(abs(data.I(i, :) - data.Ir(i, 2)));

%correct for offset voltage at I = 0;
[~, zeroIdx] = min(abs(data.I(i, :) - 0));

[f, gof] = fit(-(data.V(i, 1:(endIdx - 2))' - data.V(idx, midpt)), -data.I(i, 1:(endIdx-2))', ft,  'StartPoint', [data.RnTc, 0.025 * tanh(1.74*sqrt(80/data.Ir(i, 1) - 1))], 'problem', data.Ir(i, 1), 'Lower', [0.1, 0.0001], 'Upper', [100, 0.06], 'Display', 'final', 'TolFun', 1e-15, 'TolX', 1e-10);

figure(10);
ylabel('Current / Area (kA/cm^2)');
xlabel('Voltage (mV)');
xlim([0, 50]);
ylim([1e-2, 1])

figure(20); 
plot((0:0.001:0.05) * 1000, f(0:0.001:0.05) / data.Area / 1000 * 1e8, '--', 'color', 0.5*[1,1,1], 'linewidth', 2, 'DisplayName', 'Fit');
ylabel('Current / Area (kA/cm^2)');
xlabel('Voltage (mV)');
xlim([0, inf]);
ylim([0, 0.6]);

figure(30);
ylabel('dV/dI (\Omega)');
xlabel('Voltage (mV)');
xlim([0, inf]);
ylim([0, inf]);    

%% S7: Subgap Features in dV/dI
% 0 degree
figure(1);
hold off;
list = 1:1:17;
colors = jet(length(list));
for i = 1:length(list)
    plot(-(D0A376.V(i, :) - D0A376.V(i, 151)) * 1000, -D0A376.I(i, :) / D0A376.Area * 1e6, '.', 'color', colors(i, :));
    hold all;
end

plot([6, 6], [14, 0]);
plot([8, 8], [14, 0]);

xlim([0, 20]);
ylim([2, 14]);
ylabel('I / A (\muA / \mum^2)');
xlabel('V (mV)');


figure(2);
hold off;
list = 1:17;
colors = jet(length(list));
for i = 1:length(list)
    PlotWithBreaks(-(D0A376.V(i, 2:end) - D0A376.V(i, 151)) * 1000, D0A376.dV(i, 2:end) ./ D0A376.dI(i, 2:end) + i*0.5, 0.2, colors(i, :));
    hold all;
end
%legend for 5 ohms.
plot(18.5*[1,1], [11,13]);
ylim([2,14]);
xlim([10, 19]);
xlabel('V (mV)');
ylabel('Offset dV/dI');

%% S8: Angle dependence of various junction parameters. 

varList = whos('D*');
idx = [2	1	26	5	25	24	6	7	8	3	9	23	4	11	12	10	13	17	14	15	22	16	19	18	21	20];
varList = varList(idx);


figure(1); clf; figure(2); clf; figure(3); clf; figure(4); clf; figure(5); clf; figure(6); clf;



T = 60;
runSum = 0;
for i = 1:length(varList)
	data = eval(varList(i).name);
	midpt = ceil(size(data.V, 2) / 2);
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	[~, idx10] = min(abs(data.Ic2(:, 1) - 30));

	if size(data.I, 2) == 1
		[~, idx2] = min(abs(data.I - data.Ic2(idx, 2)));
		idx2 = idx2 + 2;

	else
		[~, idx2] = min(abs(data.I(idx, :) - data.Ic2(idx, 2)));
		idx2 = idx2 + 2;

	end

	if data.Angle <= 90
		figure(1); 
		plot(data.Angle, data.Ic2(idx, 2) / data.Area / 1000 * 1e8, 'o', 'color', [63,232,163] / 255, 'MarkerFaceColor', [63,232,163] / 255); hold all;

		figure(2); 
		plot(data.Angle, (data.V(idx, idx2) - data.V(idx, midpt)) * 1000, 'o', 'color', [255,139,0] / 255, 'MarkerFaceColor', [255,139,0] / 255); hold all;

		figure(3); 
		try
		plot(data.Angle, abs(data.Ic2(idx, 2)) * data.Rn60 * 1000, 'o', 'color', [1,0,0], 'MarkerFaceColor', [1,0,0]); hold all;
		end

		figure(4); 
		plot(data.Angle, data.Tc2, 'o', 'color', 0*[1,1,1], 'MarkerFaceColor', 0*[1,1,1]); hold all;

		figure(5);
		plot(data.Angle, 1000 / (data.RnTc2 * data.Area), 'o', 'color', 0.5*[1,1,1], 'MarkerFaceColor', 0.5*[1,1,1]); hold all;
		runSum = runSum + 1000 / (data.RnTc2 * data.Area);
	else
		figure(1); 
		plot(data.Angle - 90, data.Ic2(idx, 2) / data.Area / 1000 * 1e8, 'o', 'color', [63,232,163] / 255); hold all;

		figure(2); 
		plot(data.Angle - 90, (data.V(idx, idx2) - data.V(idx, midpt)) * 1000, 'o', 'color', [255,139,0] / 255); hold all;
				
		figure(3); 
		plot(data.Angle - 90, abs(data.Ic2(idx, 2)) * data.Rn60 * 1000, 'o', 'color', [1,0,0]); hold all;
				
		figure(4); 
		plot(data.Angle - 90, data.Tc2, 'o', 'color', 0*[1,1,1]); hold all;
		
		figure(5);
		plot(data.Angle - 90, 1000 / (data.RnTc2 * data.Area), 'o', 'color', 0.5*[1,1,1]); hold all;
		runSum = runSum + 1000 / (data.RnTc2 * data.Area);


	end
	figure(2);
	 xlim([0, 90]); ylim([0, 25]);
	
end
x = 0:90;
figure(1); 
plot(x, 0.8 * abs(cos(2*x * pi/180)), 'color', 0.35*[1,1,1])
xlabel('\theta mod 90');
ylabel('Ic / A (A/um^2)')
xticks([0, 0.25, 0.5, 0.75, 1] * 90);

figure(2); 
plot(x, 17 * abs(cos(2*x * pi/180)), 'color', 0.35*[1,1,1])
xlabel('\theta mod 90');
ylabel('V(Ic) (mV)')
xticks([0, 0.25, 0.5, 0.75, 1] * 90);

figure(3); 
plot(x, 17 * abs(cos(2*x * pi/180)), 'color', 0.35*[1,1,1])
xlabel('\theta mod 90');
ylabel('I_CR_N^{60K} (mV)')
xticks([0, 0.25, 0.5, 0.75, 1] * 90);

figure(4); 
yline(84);
xlabel('\theta mod 90');
ylabel('T_C (K)')
xticks([0, 0.25, 0.5, 0.75, 1] * 90);
ylim([0, 100]);

figure(5); 

xlabel('\theta mod 90');
ylabel('R_N^{-1} / A (mS / \mum^2)')
xticks([0, 0.25, 0.5, 0.75, 1] * 90);


%% S9: Plot angle-dependent IVCs.
colors = [46, 114, 255; 90, 224, 255; 255, 142, 136; 255, 39, 39; 148, 60, 255] / 255;

figure(1); clf; %angle dependent IV plots.
T = 7;
	data = D14A255;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot((data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.I(idx, :) / 255 * 1e5, 'color', colors(1, :)); hold all;
	disp(mean(data.TSample(idx, :)))

	data = D119A736;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot((data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.I(idx, :) / 736 * 1e5,  'color', colors(2, :)); hold all;
	disp(mean(data.TSample(idx, :)))
	
	data = D430A111;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot((data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.I(idx, :) / 111 * 1e5,  'color', colors(3, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D449A123;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot((data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.I(idx, :) / 123 * 1e5,  'color', colors(4, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D70A293;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot((data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.I(idx, :) / 293  * 1e5,  'color', colors(5, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

xlabel('V (mV)');
ylabel('J (kA/cm^2)')
xlim([-40, 40]);
ylim([-2, 2]);
legend('14', '29', '43.0', '44.9', '70')

figure(2); clf; %semilog plot of IV.
T = 7;
	data = D14A255;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, -data.I(idx, :) / 255 * 1e5, '.', 'color', colors(1, :)); hold all;
	disp(mean(data.TSample(idx, :)))

	data = D119A736;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, -data.I(idx, :) / 736 * 1e5, '.', 'color', colors(2, :)); hold all;
	disp(mean(data.TSample(idx, :)))
	
	data = D430A111;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, -data.I(idx, :) / 111 * 1e5, '.', 'color', colors(3, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D449A123;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000,-data.I(idx, :) / 123 * 1e5, '.', 'color', colors(4, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D70A293;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy( -(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, -data.I(idx, :) / 293 * 1e5, '.', 'color', colors(5, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

xlabel('V (mV)');
ylabel('J (kA/cm^2)')
xlim([0, 40]);
ylim([1e-2, 3]);

figure(3); clf; %angle dependent dI/dV plots.
T = 7;
	data = D14A255;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2, '.', 'color', colors(1, :)); hold all;
	disp(mean(data.TSample(idx, :)))

	data = D119A736;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,  '.',  'color', colors(2, :)); hold all;
	disp(mean(data.TSample(idx, :)))
	
	data = D430A111;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,  '.',  'color', colors(3, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D449A123;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,  '.',  'color', colors(4, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D70A293;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	plot(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,   '.', 'color', colors(5, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

xlabel('V (mV)');
ylabel('dV/dI / Rn')
xlim([0, 40]);
ylim([0,14]);
% legend('14', '29', '43.0', '44.9', '70')

figure(30); clf; %angle dependent dI/dV plots.
T = 7;
	data = D14A255;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2, '.', 'color', colors(1, :)); hold all;
	disp(mean(data.TSample(idx, :)))

	data = D119A736;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,  '.',  'color', colors(2, :)); hold all;
	disp(mean(data.TSample(idx, :)))
	
	data = D430A111;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,  '.',  'color', colors(3, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D449A123;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,  '.',  'color', colors(4, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

	data = D70A293;
	[~, idx] = min(abs(data.Ic2(:, 1) - T));
	semilogy(-(data.V(idx, :) - data.V(idx, ceil(end/2))) * 1000, data.dV(idx, :) ./ data.dI(idx, :) / data.RnTc2,   '.', 'color', colors(5, :)); hold all;
	disp(mean(data.TSample(idx, :))) 

xlabel('V (mV)');
ylabel('dV/dI / Rn')
xlim([0, 40]);
ylim([0.1,20]);
legend('14', '29', '43.0', '44.9', '70')


figure(4); clf; %Ic/Ir
	data = D14A255;
	plot(data.Ic2(:, 1), -4 * data.Ic2(:, 2) ./ (pi * data.Ir(:, 2)), '.', 'color', colors(1, :)); hold all;

	data = D119A736;
	plot(data.Ic2(:, 1), -4 * data.Ic2(:, 2) ./ (pi * data.Ir(:, 2)), '.', 'color', colors(2, :));
	
	data = D430A111;
	plot(data.Ic2(:, 1), -4 * data.Ic2(:, 2) ./ (pi * data.Ir(:, 2)), '.', 'color', colors(3, :));

	data = D449A123;
	plot(data.Ic2(:, 1), -4 * data.Ic2(:, 2) ./ (pi * data.Ir(:, 2)), '.', 'color', colors(4, :));

	data = D70A293;
	plot(data.Ic2(:, 1), -4 * data.Ic2(:, 2) ./ (pi * data.Ir(:, 2)), '.', 'color', colors(5, :));

	yline(4/pi);

	ylabel('4 Ic / pi Ir');
	xlabel('T (K)');
	xlim([0, 80]);

figure(5); clf; %Vc(T)
	data = D14A255;
	plot(data.Ic2(:, 1), data.Ic2(:, 3) * 1000, '.', 'color', colors(1, :)); hold all;

	data = D119A736;
	plot(data.Ic2(:, 1), data.Ic2(:, 3) * 1000, '.', 'color', colors(2, :));
	
	data = D430A111;
	plot(data.Ic2(:, 1), data.Ic2(:, 3) * 1000, '.', 'color', colors(3, :));

	data = D449A123;
	plot(data.Ic2(:, 1), data.Ic2(:, 3) * 1000, '.', 'color', colors(4, :));

	data = D70A293;
	plot(data.Ic2(:, 1), data.Ic2(:, 3) * 1000, '.', 'color', colors(5, :));

	ylabel('Vc (mV)');
	xlabel('T (K)');
	xlim([0, 80]);
	ylim([0, 35]);

figure(6); clf;
	data = D14A255;
	semilogy(data.Ic2(:, 1), (data.Ic2(:, 3) * 1000), '.', 'color', colors(1, :)); hold all;

	data = D119A736;
	semilogy(data.Ic2(:, 1), (data.Ic2(:, 3) * 1000), '.', 'color', colors(2, :));
	
	data = D430A111;
	semilogy(data.Ic2(:, 1), (data.Ic2(:, 3) * 1000), '.', 'color', colors(3, :));

	data = D449A123;
	semilogy(data.Ic2(:, 1), (data.Ic2(:, 3) * 1000), '.', 'color', colors(4, :));

	data = D70A293;
	semilogy(data.Ic2(:, 1), (data.Ic2(:, 3) * 1000), '.', 'color', colors(5, :));

	ylabel('Vc (mV)');
	xlabel('T (K)');
	xlim([0, 80]);
  	ylim([1e-2, 5e1]);

%% S10: IV curves for all junctions.
figure(1); clf;
varList = whos('D*');
idx = [2	1	26	5	25	24	6	7	8	3	9	23	4	11	12	10	13	17	14	15	22	16	19	18	21	20];
varList = varList(idx);

T = [10, 30, 50, 70, 80];

colors = [21, 89, 209;
        221, 199, 255;          
		  242, 22, 112];

colors = interp1([0; 0.5; 1], colors, linspace(0, 1, length(T))) / 255;

t = tiledlayout(7,4,'TileSpacing','Compact','Padding','Compact');

for i = 1:length(varList)
    varList(i).name
    eval(strcat('data = ', varList(i).name, ';'));

	nexttile
		for k = 1:length(T)
			[~, idx] = min(abs(data.Ic2(:, 1) - T(k)));
			if size(data.I, 2) > 1
				plot(data.I(idx, :) / data.Area / 1000 * 1e8, data.V(idx, :) * 1000, 'color', colors(k, :), 'linewidth', 2); hold all;
			else
				plot(data.I / data.Area / 1000 * 1e8, data.V(idx, :) * 1000, 'color', colors(k, :), 'linewidth', 2); hold all;
			end
		end
		ylabel('V (mV)');
		xlabel('I/A (kA/cm^2)')
		text(0.03, 0.9, sprintf('%.1f^o', data.Angle), 'Units', 'normalized')
    % pause(1)
end


%% S11: Color plots for all devices,
varList = whos('D*');
%Order runs by distance from 45*. Use abs(45 - mod(angle, 90)).
idx = [2	1	26	5	25	24	6	7	8	3	9	23	4	11	12	10	13	17	14	15	22	16	19	18	21	20];
varList = varList(idx);

t = tiledlayout(6, 4,'TileSpacing','Compact','Padding','Compact');

figure(3);
clf
a = -8; b = 8;
xMax = [20, 20, 30, 60, 60, 50, 50, 25, 30, 40, 15, 30, 30, 15, 15, 10, 10, 15, 10, 10, 10, 10, 10, 10, 10, 10];
j = 1;
for i = 1:length(varList)
%for i = 1
    nexttile;
    j = j+1;
    varList(i).name
    eval(strcat('data = ', varList(i).name, ';'));
    if size(data.I, 2) == 1
        I = repmat(data.I', size(data.V, 1), 1);
    else
        I = data.I;
    end

    if ~isfield(data, 'dVdI') && ~isfield(data, 'CdVdI')
        data.dVdI = data.dV ./ data.dI;
    end

    try
        pcolor(I(:, 2:end) * data.RnTc2 * 1e3, data.TSample(:, 2:end), data.CdVdI / data.RnTc2);
        hold all;

    

    catch
        pcolor(I(:, 1:end) * data.RnTc2 * 1e3, data.TSample(:, 1:end), data.dVdI / data.RnTc2);
        hold all;

    end
        shading interp;
        ylim([4, 90]);
        xlim([-xMax(i), xMax(i)]);
        caxis([0, 4]);
        map = GenerateColorScaleJunctions();
        colormap(map);

        hold all
    ylabel('Temperature (K)');
    xlabel('I Rn (mV)');

    title(strcat('\theta = ', num2str(data.Angle), '^\circ, \phi = ', num2str(mod(data.Angle, 90)), '^\circ, \alpha = ', num2str(data.alpha), '^\circ'));
end


%% S13: Ic Repeatability
colors = winter(7);
figure(1); clf;
plot(IcRepeatA.a(1, :), IcRepeatA.a(2, :) / 446 / 1000 * 1e8, '.'); hold all;%, 'color', colors(1, :)); hold all;
 plot(IcRepeatA.b(1, :), IcRepeatA.b(2, :) / 446 / 1000 * 1e8, '.'); %, 'color', colors(3, :));
 plot(IcRepeatA.c(1, :), IcRepeatA.c(2, :) / 446 / 1000 * 1e8, '.') %, 'color', colors(5, :));

%main dataset.
plot(D39A446.Ic2(:, 1), D39A446.Ic2(:, 2) / 446 / 1000 * 1e8, '.', 'color', [0,0,0])

xlabel('Temperature (K)')
ylabel('Jc (kA/cm^2)');

% pause(1);
figure(2); clf;
plot(IcRepeatB.a(1, :), IcRepeatB.a(2, :) / 417 / 1000 * 1e8, '.'); hold all;% , 'color', colors(3, :)); hold all;
plot(IcRepeatB.b(1, :), IcRepeatB.b(2, :) / 417 / 1000 * 1e8, '.'); %, 'color', colors(4, :));
plot(IcRepeatB.c(1, :), IcRepeatB.c(2, :) / 417 / 1000 * 1e8, '.'); %, 'color', colors(5, :));
plot(D31A417.Ic2(:, 1), D31A417.Ic2(:, 2) / 417 / 1000 * 1e8, '.', 'color', [0,0,0])

xlabel('Temperature (K)')
ylabel('Jc (kA/cm^2)');
